{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import beta\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize']=(15,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "a =2;b=2    # 先验分布中的伪计数\n",
    "N1 = 3;N0 = 17   # 观察到的实验结果\n",
    "M = 10     # 需要预测的实验次数 ，请参考文中的字母表达"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'(b)')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAAFBCAYAAAA/uO2MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGelJREFUeJzt3X+MZedZH/Dvw26dQFKoIVto1t5kI5YGU7W4GRzatOFHfjmi9fJHUBx+yJRIK0rc0qZV5ZSq2TWiClBVRcKCWOCCKI0JbgIrFDCJE0BIOOxsEiWxg5vNJtiDQx1wIG0T4mzy9I+5huuZ2Z2769k59939fKTR3nPOe+4+99g7z3znvOec6u4AAAAwni+augAAAAAujEAHAAAwKIEOAABgUAIdAADAoAQ6AACAQQl0AAAAg1oo0FXV9VX1QFWdqqpbttj+/VX1gap6X1X9blVdM7ftdbP9Hqiql+1k8QAwNT0SgCnVds+hq6o9Sf5XkpckWUtyIsmruvv+uTFf2t2fmr2+IckPdPf1s6b1piTXJXlmknck+Zru/vzF+DAAsJv0SACmtsgZuuuSnOru0939WJI7kxyeH/B4o5p5WpLHU+LhJHd292e7+6NJTs3eDwAuBXokAJPau8CY/UkemlteS/L8jYOq6jVJXpvkiiTfOrfvvRv23X9BlQLA8tEjAZjUIoGutli3aZ5md9+W5Laq+s4k/yHJTYvuW1VHkhxJkqc97WnPe+5zn7tAWQCM7uTJk3/S3fumruNJ0CO3cfLkhe/7vOftXB0AIzmf/rhIoFtLcvXc8lVJHj7H+DuT/NT57Nvdtye5PUlWVlZ6dXV1gbIAGF1V/eHUNTxJeuQ2aqvYuqDBPirAjjmf/rjINXQnkhyqqoNVdUWSG5Mc3/AXHppb/LYkH569Pp7kxqp6SlUdTHIoye8vWhwALDk9EoBJbXuGrrvPVNXNSe5OsifJHd19X1XdmmS1u48nubmqXpzkc0k+mfWpJJmNe3OS+5OcSfIad+8C4FKhRwIwtW0fW7DbRpxOAsCFqaqT3b0ydR2jGLFHPpkplzv5I8rRo9PsC3Ahzqc/LnINHQDA0I4du/B9BTpgmS1yDR0AAABLSKADAAAYlCmXO8C8fAAAYAoC3Q4wLx8AAJiCKZcAAACDEugAAAAGJdABAAAMSqADAAAYlEAHAAAwKIEOAABgUAIdAADAoAQ6AACAQXmw+CXmyTyo3EPOAQBgLALdJebYsQvfV6ADAICxmHIJAAAwKIEOAABgUAIdAADAoAQ6AACAQQl0AAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACDWijQVdX1VfVAVZ2qqlu22P7aqrq/qt5fVfdU1bPmtn2+qt43+zq+k8UDwJT0RwCmtne7AVW1J8ltSV6SZC3Jiao63t33zw17b5KV7v50Vf3zJD+W5JWzbZ/p7q/f4boBYFL6IwDLYJEzdNclOdXdp7v7sSR3Jjk8P6C739Xdn54t3pvkqp0tEwCWjv4IwOQWCXT7kzw0t7w2W3c2r07y63PLT62q1aq6t6q+/QJqBIBlpD8CMLltp1wmqS3W9ZYDq747yUqSb5pbfaC7H66q5yR5Z1V9oLs/smG/I0mOJMmBAwcWKhwAJnbR++NsXz0SgLNa5AzdWpKr55avSvLwxkFV9eIkP5Tkhu7+7OPru/vh2Z+nk/xWkms37tvdt3f3Snev7Nu377w+AABM5KL3x9l2PRKAs1ok0J1IcqiqDlbVFUluTPKEu3FV1bVJ3pj1ZvXI3Porq+ops9fPSPKCJPMXiwPAqPRHACa37ZTL7j5TVTcnuTvJniR3dPd9VXVrktXuPp7kx5M8PckvV1WSPNjdNyT52iRvrKovZD08vmHD3b8AYEj6IwDLoLq3nO4/mZWVlV5dXZ26jPNSW11FsaCdPvzLVAvAdqrqZHevTF3HKPTI8esAWMT59MeFHiwOAADA8hHoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACDEugAAAAGJdABAAAMSqADAAAYlEAHAAAwKIEOAABgUAIdAADAoAQ6AACAQQl0AAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADCohQJdVV1fVQ9U1amqumWL7a+tqvur6v1VdU9VPWtu201V9eHZ1007WTwATEl/BGBq2wa6qtqT5LYkL09yTZJXVdU1G4a9N8lKd//dJHcl+bHZvl+e5PVJnp/kuiSvr6ord658AJiG/gjAMljkDN11SU519+nufizJnUkOzw/o7nd196dni/cmuWr2+mVJ3t7dj3b3J5O8Pcn1O1M6AExKfwRgcosEuv1JHppbXputO5tXJ/n1C9wXAEahPwIwub0LjKkt1vWWA6u+O8lKkm86n32r6kiSI0ly4MCBBUoCgMld9P4421ePBOCsFjlDt5bk6rnlq5I8vHFQVb04yQ8luaG7P3s++3b37d290t0r+/btW7R2AJjSRe+PiR4JwLktEuhOJDlUVQer6ookNyY5Pj+gqq5N8sasN6tH5jbdneSlVXXl7GLvl87WAcDo9EcAJrftlMvuPlNVN2e90exJckd331dVtyZZ7e7jSX48ydOT/HJVJcmD3X1Ddz9aVT+c9aaXJLd296MX5ZMAwC7SHwFYBtW95ZT9yaysrPTq6urUZZyX2upKiAXt9OFfploAtlNVJ7t7Zeo6RqFHjl8HwCLOpz8u9GBxAAAAlo9ABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACDEugAAAAGJdABAAAMSqADAAAYlEAHAAAwKIEOAABgUHunLoBL19Gj0+wLAACXC4GOi+bYsQvfV6ADAIDtmXIJAAAwKIEOAABgUKZcAgDsIteYAztJoAMA2EWuMQd2kimXAAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxqoUBXVddX1QNVdaqqbtli+wur6j1VdaaqXrFh2+er6n2zr+M7VTgATE1/BGBqe7cbUFV7ktyW5CVJ1pKcqKrj3X3/3LAHk3xvkn+7xVt8pru/fgdqBYCloT8CsAy2DXRJrktyqrtPJ0lV3ZnkcJK/bFjd/bHZti9chBoBYBnpjwBMbpEpl/uTPDS3vDZbt6inVtVqVd1bVd9+XtUBwPLSHwGY3CJn6GqLdX0ef8eB7n64qp6T5J1V9YHu/sgT/oKqI0mOJMmBAwfO460BYDIXvT8meiQA57bIGbq1JFfPLV+V5OFF/4Lufnj25+kkv5Xk2i3G3N7dK929sm/fvkXfGgCmdNH742y7HgnAWS1yhu5EkkNVdTDJHyW5Mcl3LvLmVXVlkk9392er6hlJXpDkxy602I2OHp1mXwDIEvdHAC4f2wa67j5TVTcnuTvJniR3dPd9VXVrktXuPl5V35DkrUmuTPJPq+pYd39dkq9N8sbZxeBflOQNG+7+9aQcO3bh+wp0ADwZy9wfAbh8VPf5TPe/+FZWVnp1dXWhsbXV1QsL2smPvSx1JGoBxlJVJ7t7Zeo6RnE+PXJZLEsvWJY6kuWqBVhO59MfF5lyCQAMxCUJAJcPgQ4ALjEuSQC4fCxyl0sAAACWkEAHAAAwKIEOAABgUAIdAADAoAQ6AACAQQl0AAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACDEugAAAAGJdABAAAMSqADAAAYlEAHAAAwKIEOAABgUAIdAADAoAQ6AACAQQl0AAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMKiFAl1VXV9VD1TVqaq6ZYvtL6yq91TVmap6xYZtN1XVh2dfN+1U4QAwNf0RgKltG+iqak+S25K8PMk1SV5VVddsGPZgku9N8j827PvlSV6f5PlJrkvy+qq68smXDQDT0h8BWAaLnKG7Lsmp7j7d3Y8luTPJ4fkB3f2x7n5/ki9s2PdlSd7e3Y929yeTvD3J9TtQNwBMTX8EYHKLBLr9SR6aW16brVvEk9kXAJaZ/gjA5BYJdLXFul7w/Rfat6qOVNVqVa1+4hOfWPCtAWBSF70/JnokAOe2SKBbS3L13PJVSR5e8P0X2re7b+/ule5e2bdv34JvDQCTuuj9MdEjATi3RQLdiSSHqupgVV2R5MYkxxd8/7uTvLSqrpxd7P3S2ToAGJ3+CMDktg103X0myc1ZbzQfSvLm7r6vqm6tqhuSpKq+oarWknxHkjdW1X2zfR9N8sNZb3onktw6WwcAQ9MfAVgG1b3odP/dsbKy0qurqwuNra2uQFjQTn7sZakjUQswlqo62d0rU9cxikV75DJ9/12WWpaljmS5agGW0/n0x4UeLA4AAMDyEegAAAAGJdABAAAMSqADAAAYlEAHAAAwKIEOAABgUAIdAADAoAQ6AACAQQl0AAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIPaO3UBcLEdPTrNvgAAcLEJdFzyjh278H0FOgAAlpkplwAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACDEugAAAAGtXfqAgAAmMbRo9PsC+wcgQ4A4DJ17NiF7yvQwXIw5RIAAGBQAh0AAMCgBDoAAIBBCXQAAACDEugAAAAGJdABAAAMaqFAV1XXV9UDVXWqqm7ZYvtTquqXZtvfXVXPnq1/dlV9pqreN/v66Z0tHwCmoz8CMLVtn0NXVXuS3JbkJUnWkpyoquPdff/csFcn+WR3f3VV3ZjkR5O8crbtI9399TtcNwBMSn8EYBkscobuuiSnuvt0dz+W5M4khzeMOZzk52ev70ryoqqqnSsTAJaO/gjA5BYJdPuTPDS3vDZbt+WY7j6T5M+TfMVs28Gqem9V/XZV/eMnWS8ALAv9EYDJbTvlMslWv0nsBcd8PMmB7v7Tqnpekl+pqq/r7k89YeeqI0mOJMmBAwcWKAkAJnfR+2OiRwJwboucoVtLcvXc8lVJHj7bmKram+TLkjza3Z/t7j9Nku4+meQjSb5m41/Q3bd390p3r+zbt+/8PwUA7L6L3h9n2/VIAM5qkUB3IsmhqjpYVVckuTHJ8Q1jjie5afb6FUne2d1dVftmF42nqp6T5FCS0ztTOgBMSn8EYHLbTrns7jNVdXOSu5PsSXJHd99XVbcmWe3u40l+NskvVNWpJI9mvaklyQuT3FpVZ5J8Psn3d/ejF+ODAMBu0h8BWAbVvXG6/7RWVlZ6dXV1obFP5j5hO/mxl6WORC3LXAewWVWd7O6VqesYxaI9cpm+7y1LLctSR6IWYHvn0x8XerA4AAAAy0egAwAAGJRABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6AAAAAYl0AEAAAxKoAMAABiUQAcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACDEugAAAAGJdABAAAMSqADAAAY1N6pC4DLydGj0+wLAMClSaCDXXTs2IXvK9ABALCRKZcAAACDEugAAAAGJdABAAAMSqADAAAYlEAHAAAwKHe5BABgUh7rAxdOoAMAYFIe6wMXzpRLAACAQQl0AAAAgxLoAAAABiXQAQAADEqgAwAAGJRABwAAMCiBDgAAYFCeQweXIQ9wBQC4NAh0cBnyAFcAgEvDQlMuq+r6qnqgqk5V1S1bbH9KVf3SbPu7q+rZc9teN1v/QFW9bOdKB4Dp6ZEATGnbM3RVtSfJbUlekmQtyYmqOt7d988Ne3WST3b3V1fVjUl+NMkrq+qaJDcm+bokz0zyjqr6mu7+/E5/EADYbXokXHpclsBoFplyeV2SU919Okmq6s4kh5PMN6vDSY7OXt+V5Cerqmbr7+zuzyb5aFWdmr3f7+1M+QAwKT0SLjEuS2A0iwS6/UkemlteS/L8s43p7jNV9edJvmK2/t4N++6/4GqBS47fhDI4PRK4KPRHFrVIoKst1vWCYxbZN1V1JMmR2eL/raoHFqhrO89I8idn21hbVTaBXa5jiGOS7GotjslmwxyTJ/Nb1PN0zmNyGduJ4/KsnShkQiP2yGH+jfu+t5mfGzZzTDbbxf6Y6JFb2dX+uEigW0ty9dzyVUkePsuYtaram+TLkjy64L7p7tuT3L5o0YuoqtXuXtnJ9xydY7KZY7KZY7KZY7I1xyXJgD3Sf7fNHJOtOS6bOSabOSab7fYxWeQulyeSHKqqg1V1RdYv4D6+YczxJDfNXr8iyTu7u2frb5zd4etgkkNJfn9nSgeAyemRAExq2zN0s/n+Nye5O8meJHd0931VdWuS1e4+nuRnk/zC7ILuR7Pe0DIb9+asXxx+Jslr3L0LgEuFHgnA1BZ6sHh3vy3J2zas+49zr/8iyXecZd8fSfIjT6LGC7WjUzgvEY7JZo7JZo7JZo7J1hyXDNkj/XfbzDHZmuOymWOymWOy2a4ek1qf9QEAAMBoFrmGDgAAgCV0yQW6qrq+qh6oqlNVdcvU9Uytqq6uqndV1Yeq6r6q+sGpa1oWVbWnqt5bVb82dS3Loqr+RlXdVVV/MPt/5h9MXdPUqupfz/7tfLCq3lRVT526pt1WVXdU1SNV9cG5dV9eVW+vqg/P/rxyyhpZjB75RHrk2emRT6Q/bqY/rluGHnlJBbqq2pPktiQvT3JNkldV1TXTVjW5M0n+TXd/bZJvTPIax+Qv/WCSD01dxJL5iSS/0d3PTfL3cpkfn6ran+RfJlnp7r+T9Zte3DhtVZP4uSTXb1h3S5J7uvtQkntmyywxPXJLeuTZ6ZFPpD/O0R+f4OcycY+8pAJdkuuSnOru0939WJI7kxyeuKZJdffHu/s9s9f/J+vfgPZPW9X0quqqJN+W5GemrmVZVNWXJnlh1u/Il+5+rLv/bNqqlsLeJF88e37Yl2SL54Rd6rr7d7J+d8Z5h5P8/Oz1zyf59l0tiguhR26gR25Nj3wi/fGsLvv+mCxHj7zUAt3+JA/NLa/FN+a/VFXPTnJtkndPW8lS+K9J/l2SL0xdyBJ5TpJPJPlvs2k2P1NVT5u6qCl19x8l+c9JHkzy8SR/3t2/OW1VS+Mru/vjyfoPxUn+5sT1sD098hz0yCfQI59If9xAf9zWrvbISy3Q1Rbr3MYzSVU9Pcn/TPKvuvtTU9czpar6J0ke6e6TU9eyZPYm+ftJfqq7r03y/3KZT6ObzXk/nORgkmcmeVpVffe0VcEF0yPPQo/8K3rklvTHDfTH5XKpBbq1JFfPLV+Vy/T077yq+mtZb1S/2N1vmbqeJfCCJDdU1ceyPuXoW6vqv09b0lJYS7LW3Y//dvqurDewy9mLk3y0uz/R3Z9L8pYk/3DimpbF/66qv5Uksz8fmbgetqdHbkGP3ESP3Ex/3Ex/PLdd7ZGXWqA7keRQVR2sqiuyfnHm8YlrmlRVVdbnfH+ou//L1PUsg+5+XXdf1d3Pzvr/I+/s7sv+t0rd/cdJHqqqvz1b9aIk909Y0jJ4MMk3VtWXzP4tvSiX+YXwc44nuWn2+qYkvzphLSxGj9xAj9xMj9xMf9yS/nhuu9oj917MN99t3X2mqm5OcnfW77ZzR3ffN3FZU3tBku9J8oGqet9s3b/v7rdNWBPL618k+cXZD3unk/yzieuZVHe/u6ruSvKerN8N771Jbp+2qt1XVW9K8s1JnlFVa0len+QNSd5cVa/OemP/jukqZBF65Jb0SBalP87RH//KMvTI6jZ9HgAAYESX2pRLAACAy4ZABwAAMCiBDgAAYFACHQAAwKAEOgAAgEEJdAAAAIMS6GACVfXFVfXbVbXnLNuvqKrfqapL6lmRAHAuc/3xRVX1a2cZ846qunK3a4NlJdDBNL4vyVu6+/Nbbezux5Lck+SVu1oVAEzr+5K8JcmW/XHmF5L8wO6UA8tPoINpfFeSX62qp1fVPVX1nqr6QFUdnhvzK7NxAHC5+K4kvzp7/aVV9daqur+qfrqqHv+59XiSV01THiyf6u6pa4DLSlVdkeTB7v6q2ZTKL+nuT1XVM5Lcm+RQd/dsOuYfd/e+SQsGgF2woT9+c5LfSHJNkj+cvX5jd981G/vhJN/Y3X86Vb2wLJyhg933jCR/NntdSf5TVb0/yTuS7E/ylUkym475WFX99UmqBIDdNd8fk+T3u/v0rB++Kck/mtv2SJJn7mZxsKzccAF232eSPHX2+ruS7EvyvO7+XFV9bG5bkjwlyV/sbnkAMIn5/pgkG6eRzS8/dTYeLnvO0MEu6+5PJtlTVU9N8mVJHpmFuW9J8qzHx1XVVyT5RHd/bqJSAWDXbOiPSXJdVR2cXTv3yiS/myRVVUm+KsnHJikUloxAB9P4zaxPHfnFJCtVtZr1s3V/MDfmW5K8bYLaAGAqj/fHJPm9JG9I8sEkH03y1tn65yW5t7vP7H55sHzcFAUmUFXXJnltd3/POca8JcnruvuB3asMAKazYH/8iSTHu/ue3asMlpczdDCB7n5vkned68HiSX5FmAPgcrJdf5z5oDAHf8UZOgAAgEE5QwcAADAogQ4AAGBQAh0AAMCgBDoAAIBBCXQAAACD+v8KlOkWUQLucQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 使用贝叶斯模型平均进行预测\n",
    "ax1 = plt.subplot(121)\n",
    "ax2 = plt.subplot(122)\n",
    "a = a+N1;b=b+N0\n",
    "result = np.empty(M+1)\n",
    "for x in range(M+1):\n",
    "    result[x]=np.math.factorial(M)/(np.math.factorial(x)*np.math.factorial(M-x))* \\\n",
    "    (beta(x+a,M-x+b)/beta(a,b))\n",
    "ax1.vlines(np.arange(M+1),0,result,colors='blue',linewidth=15)\n",
    "ax1.set_ylim((0,0.3))\n",
    "ax1.set_xlabel('(a)')\n",
    "\n",
    "# 基于最大后验估计的点估计\n",
    "plugin = (a-1)/(a+b-2)\n",
    "result_2 = np.empty(M+1)\n",
    "for x in range(M+1):\n",
    "    result_2[x]= np.math.factorial(M)/(np.math.factorial(x)*np.math.factorial(M-x))*np.power(plugin,x)*np.power(1-plugin,M-x)\n",
    "ax2.vlines(np.arange(M+1),0,result_2,colors='blue',linewidth=15)\n",
    "ax2.set_ylim((0,0.3))\n",
    "ax2.set_xlabel('(b)')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
